{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "134e7f9d",
   "metadata": {},
   "source": [
    "# Example 5: Special functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2571d531",
   "metadata": {},
   "source": [
    "Let's construct a dataset which contains special functions $f(x,y)={\\rm exp}(J_0(20x)+y^2)$, where $J_0(x)$ is the Bessel function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2075ef56",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cuda\n",
      "checkpoint directory created: ./model\n",
      "saving model version 0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "| train_loss: 5.15e-01 | test_loss: 5.86e-01 | reg: 5.84e+00 | : 100%|█| 20/20 [00:03<00:00,  5.89it\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "saving model version 0.1\n"
     ]
    }
   ],
   "source": [
    "from kan import *\n",
    "\n",
    "device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n",
    "print(device)\n",
    "\n",
    "# create a KAN: 2D inputs, 1D output, and 5 hidden neurons. cubic spline (k=3), 5 grid intervals (grid=5).\n",
    "model = KAN(width=[2,1,1], grid=3, k=3, seed=2, device=device)\n",
    "f = lambda x: torch.exp(torch.special.bessel_j0(20*x[:,[0]]) + x[:,[1]]**2)\n",
    "dataset = create_dataset(f, n_var=2, device=device)\n",
    "\n",
    "# train the model\n",
    "model.fit(dataset, opt=\"LBFGS\", steps=20);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f30c3ab",
   "metadata": {},
   "source": [
    "Plot trained KAN, the bessel function shows up in the bettom left"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3f95fcdd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwjUlEQVR4nO3deXRUZZ4+8OetqiyVjSwkgRAQEoslQFjCEkOQoLSxwW3Egw52HwF7RtGWEafnDK22IC2Ltj0SGpnpxvGo7QI94NFGEBCaPWAgBIIsWSwxJCF7KlslqeW+vz8k90c0hCy3UpXk+ZzjH52qm/qSrreeetcrpJQSREREGtK5uwAiIup7GC5ERKQ5hgsREWmO4UJERJpjuBARkeYYLkREpDmGCxERaY7hQkREmmO4EBGR5hguRESkOYYLERFpjuFCRESaY7gQEZHmGC5ERKQ5hgsREWnO4O4CiHoDKSUqKytRX1+PgIAAhIWFQQjh7rKIPBZ7LkTtsFgsSEtLg8lkQnh4OEaMGIHw8HCYTCakpaXBYrG4u0QijyR4J0qitu3duxfz58+H1WoF8EPvpUVLr8XPzw87duxAamqqW2ok8lQMF6I27N27F/PmzYOUEoqi3PR5Op0OQgjs2rWLAUN0A4YL0Y9YLBZER0ejsbGx3WBpodPpYDQaUVhYiODgYNcXSNQLcM6F6Efef/99WK3WDgULACiKAqvVig8++MDFlRH1Huy5EN1ASgmTyQSz2YzONA0hBGJiYpCXl8dVZERguBC1UlFRgfDw8G5dHxYWpmFFRL0Th8WIblBfX9+t6+vq6jSqhKh3Y7gQ3SAgIKBb1wcGBmpUCVHvxnAhukFYWBhiY2M7PW8ihEBsbCxCQ0NdVBlR78JwIbqBEALPPfdcl65dtmwZJ/OJruOEPtGPcJ8LUfex50L0I8HBwdixYweEENDp2m8iLTv0P/30UwYL0Q0YLkRtSE1Nxa5du2A0GiGE+MlwV8vPjEYjdu/ejXvuucdNlRJ5JoYL0U2kpqaisLAQGzZsQExMTKvHYmJisGHDBhQVFTFYiNrAOReiDpBS4uDBg7j77rtx4MABzJ49m5P3RO1gz4WoA4QQ6pxKcHAwg4XoFhguRESkOYYLERFpjuFCRESaY7gQEZHmGC5ERKQ5hgsREWmO4UJERJpjuBARkeYYLkREpDmGCxERaY7hQkREmmO4EBGR5hguRESkOYYLERFpjuFCRESaY7gQEZHmGC5Et2C321FUVIRLly4BAL799ltUVVVBURQ3V0bkuXibY6KbsFgs2LFjBz766CNcuHABdXV1sNls8PX1RXh4OGbOnIknn3wSM2bMgMFgcHe5RB6F4ULUhhMnTmD58uXIzs7G1KlTMW/ePMTHxyMgIAAWiwWZmZnYuXMn8vPz8eijj+K1115DeHi4u8sm8hgMF6If2bdvHxYtWoSAgACsW7cOc+fOhc1mw9atW9Hc3IygoCA89thjsNvt2Lp1K1atWoWxY8fir3/9KyIjI91dPpFHYLgQ3SA3Nxf33nsv/P39sXXrVsTFxUEIAbPZjMmTJ6OmpgYjRoxAZmYmQkJCIKXEsWPHsHDhQqSkpOCdd96Bj4+Pu/8ZRG7HCX2i65xOJ9auXYvq6mps2rRJDZb2CCGQnJyMN954A59//jn27NnTQ9USeTaGC9F1+fn52LlzJx5++GEkJyffMlhaCCHw0EMPITExEVu2bIHD4XBxpUSej0tciK5LT09HfX095s+fjytXrqChoUF9rLCwEE6nEwBgs9lw4cIFBAUFqY9HRUXh4YcfxqpVq1BSUoLo6Oger5/IkzBciK67fPky/Pz8EBMTg6eeegrHjx9XH5NSorm5GQBQXFyMn/3sZ+pjQgj88Y9/xPjx42G1WlFcXMxwoX6P4UJ0XWNjIwwGA3x8fNDc3IympqY2nyel/MljDocDRqOxVQgR9WcMF6LrIiIi0NjYCIvFgunTp8Pf3199rLGxEenp6WqIJCUlqRsnhRAYNmwYysrKoNPpEBIS4q5/ApHHYLgQXZeQkAC73Y6MjAy8/vrrrR4zm82YOnUqampqEBkZiW3btiE4OFh9XAiBF198EYMGDeKQGBG4WoxINW3aNMTExOD9999HQ0MD9Hp9q/9aCCGg0+nUn+t0Oly7dg3bt2/HvHnzMGDAADf+K4g8A8OF6LqwsDD8+te/xpkzZ7Bx48YOLylubm7G73//ezQ2NuKpp57q8BJmor6Mw2JEN1i0aBGOHDmC119/HX5+fli6dCl8fX0BAAaDAQaDQe3FSClRV1eHNWvWYOvWrXjrrbcwatQod5ZP5DF4/AvRj5SXl+PZZ5/FF198gdTUVCxfvhxjxoxBTk4OFEWBt7c3br/9dmRkZODNN9/E2bNnsXr1aixdurTV8BlRf8ZwIWpDQ0MDtmzZgo0bN6K0tBQxMTEwmUwIDAxEdXU1cnJyUFxcjISEBKxcuRKzZs2CTsdRZqIWDBeidpSUlODAgQM4fPgwzp07h4yMDMycORMzZszAPffcg+nTp8PPz8/dZRJ5HIYLUQedOnUK06ZNw6lTpzBlyhR3l0Pk0diPJ+ogvV6vLkMmovaxlRARkeYYLkREpDmGCxERaY7hQkREmmO4EBGR5hguRESkOYYLERFpjuFCRESaY7gQEZHmGC5ERKQ5hgsREWmO4UJERJpjuBARkeYYLkREpDnez4Wog6SUUBQFOp0OQgh3l0Pk0dhzIeoE3suFqGMM7i6ASAt2ux0FBQVQFMXdpXSbEALDhg2Dt7e3u0sh6jKGC/UJhYWFeOaZZ5CQkODuUrrMarXCz88PmZmZ2Lx5M2JjY91dElGXMVyoT5BSIj4+HmvWrHF3KV1y6NAhrFq1Cv/1X/8Fp9MJToVSb8dwoT6nN022OxwObNu2DS+88AIqKirw8ssvY+TIke4ui6jbODtJ5AZSSlitVqxZswZPPfUUKioqMGHCBKxbtw4+Pj7uLo+o2xguRD1MSomioiIsXrwYa9asQVNTE+bOnYvPPvsMEyZM6FU9L6KbYbgQ9SBFUXD48GHMmzcP27dvh16vx7PPPosPP/wQQ4cOdXd5RJrhnAtRD2gZBtu0aRNef/111NTUIDw8HGvXrsUvfvELeHl5QQjBiXzqMxguRC4mpURubi5+85vfYM+ePZBSYvr06UhLS8OUKVM4DEZ9EofFiFzIbrdj27ZtSE1Nxe7du+Ht7Y1nnnkGf//73xks1Kex50LkAlJKlJeXY9WqVXjvvffQ3NyM4cOH4/XXX8eDDz4Ig8HAYKE+jeFCpDFFUZCRkYFly5YhMzMTer0eDz74IN544w3ExsYyVKhfYLgQaURKCbvdjvfeew8vv/wyKisrERwcjJdeeglPP/00jEYjg4X6DYYLkQaklLBYLHj55Zfxv//7v7Db7YiPj8fGjRsxY8YMnqZM/Q7DhaibpJS4cuUKli5div3790On0+Gxxx7DH/7wBwwePJi9FeqXGC5E3SClRFZWFhYvXoxvvvkG/v7++O1vf4vnn38evr6+DBbqtxguRF0kpcThw4exePFiFBQUICIiAmlpaZg/fz70er27yyNyK4YLURdIKfHVV19h8eLFKCkpwYgRI/Duu+9i5syZ7K0QgeFC1GlSShw4cACLFi1CaWkpxowZgw8//JCHThLdgOFC1AlSSpw4cQJLlixBaWkp4uLisHXrVsTFxTFYiG7A9ZFEHSSlxKVLl7B48WIUFRXBZDLh448/ZrAQtYHhQtQBUkqUlpbiX/7lX5Cfn48hQ4bg/fffx7hx4xgsRG1guBB1QGNjI1544QWcPHkSAwYMwP/8z/9g2rRpDBaim2C4EN2C0+nExo0bsX37dnh7e2PNmjW49957GSxE7WC4ELVDSon9+/dj/fr1UBQFS5YswZNPPsnjXIhugS2E6CaklCgsLMQLL7yAuro63HHHHVi9ejW8vLzcXRqRx2O4EN2EzWbDSy+9hMuXLyMiIgIbNmxAaGgoh8OIOoDhQtQGKSW2bduGv/3tbzAYDFi5ciUmT57MYCHqIIYL0Y9IKZGfn4+VK1fCbrfjwQcfxBNPPMFgIeoEhgvRj9hsNvzud79DQUEBhg0bhrVr18LX19fdZRH1KgwXohtIKfF///d/+Oyzz+Dl5YVXX32VtyYm6gKGC9F1UkoUFBTg1Vdfhd1uxwMPPIAFCxYwWIi6gOFCdJ3T6cSaNWtgNpsRFRWF1atXw8fHx91lEfVKDBci/P/Nkh9//DH0ej3+8z//E6NGjWKvhaiLGC5EAKqrq/HKK6+gsbERd955JxYtWsRgIeoGhgv1e4qiYPPmzThz5gyCgoKwevVq+Pv7u7ssol6N4UL9mpQS2dnZ2LhxI6SU+Nd//VdMnz6dvRaibmK4UL/W1NSElStXorKyEnFxcXjhhReg1+vdXRZRr8dwoX5LSom//e1v2LNnD7y9vbFy5UpERES4uyyiPoHhQv1Sy56W1157DQ6HA//0T/+EBx54gMNhRBphuFC/5HA41D0tgwcPxsqVK3mUPpGGGC7U70gp8eWXX7ba0zJy5Ej2Wog0xHChfkVKiaKiIrz44otobGzE7NmzsXjxYgYLkcYYLtSvNDc346WXXsKlS5cQERGB9evXc08LkQswXKjPkVJCSvmTnzudTvz3f/83tm7dCoPBgFdeeQUTJ05kr4XIBRgu1Kc0NTVh27ZtqKuraxUwiqJg27ZtWLVqFRwOBx5//HEOhxG5EMOF+gyHw4GVK1di0aJFePrpp1FWVgYpJaxWKzZv3oxnn30W9fX1mD17Nl5//XWeeEzkQgZ3F0CkFSEEAgICIITAtm3bcPbsWSQlJeHSpUs4ffo0HA4HZs6ciXfffRcDBw5kr4XIhRgu1Gfo9XqsWLECISEh+P3vf4+cnBzk5OQAAHx9ffH4449j3bp1iIiIYLAQuRjDhfoULy8vPPvss7jrrruwfft25OXlISoqCvfffz8SExNhMPzwlm9rwp+ItMNwoT5BCIHz58/j1VdfbfXz2NhYCCFw4MABHDhwwE3Vdc65c+fYs6JeT0h+haM+wGazwWw2w+l0uruUbtPpdIiNjYW3t7e7SyHqMoYLERFpjsNiRB104/cwDlsRtY/7XIg6KCsrC3q9HllZWe4uhcjjMVyIiEhzDBciItIcw4WIiDTHcCEiIs0xXIiISHMMFyIi0hzDhYiINMdwISIizTFciIhIcwwXIiLSHMOFiIg0x3AhIiLNMVyIiEhzDBciItIcw4WIiDTHcCHqACklqqurAQDV1dXgDVyJ2sdwIWqHxWJBWloaTCYT5syZAykl5syZA5PJhLS0NFgsFneXSOSRhORXMKI27d27F/Pnz4fVagXQ9m2O/fz8sGPHDqSmprqlRiJPxXAhasPevXsxb948SCmhKMpNn6fT6SCEwK5duxgwRDdguBD9iMViQXR0NBobG9sNlhY6nQ5GoxGFhYUIDg52fYFEvQDnXIh+5P3334fVau1QsACAoiiwWq344IMPXFwZUe/BngvRDaSUMJlMMJvNnVoRJoRATEwM8vLy1PkYov6M4UJ0g4qKCoSHh3fr+rCwMA0rIuqdOCxGdIP6+vpuXV9XV6dRJUS9G8OF6AYBAQHduj4wMFCjSoh6N4YL0Q3CwsIQGxvb6XkTIQRiY2MRGhrqosqIeheGC9ENhBB47rnnunTtsmXLOJlPdB0n9Il+hPtciLqPPReiHwkODsaOHTsghIBO134Tadmh/+mnnzJYiG7AcCFqQ2pqKnbt2gWj0QghxE+Gu1p+ZjQasXv3btxzzz1uqpTIMzFciG4iNTUVhYWF2LBhA2JiYlo9FhMTgw0bNqCoqIjBQtQGzrkQdYCUEgcPHsTdd9+NAwcOYPbs2Zy8J2oHey5EHSCEUOdUgoODGSxEt8BwISIizTFciIhIcwwXIiLSHMOFiIg0x3AhIiLNMVyIiEhzDBciItIcw4WIiDTHcCEiIs0xXIiISHMMFyIi0hzDhYiINMdwISIizTFciIhIcwwXIiLSHMOFiIg0x3AhugW73Y6ioiJcunQJAPDtt9+iqqoKiqK4uTIiz8XbHBPdhMViwY4dO/DRRx/hwoULqKurg81mg6+vL8LDwzFz5kw8+eSTmDFjBgwGg7vLJfIoDBeiNpw4cQLLly9HdnY2pk6dinnz5iE+Ph4BAQGwWCzIzMzEzp07kZ+fj0cffRSvvfYawsPD3V02kcdguBD9yL59+7Bo0SIEBARg3bp1mDt3Lmw2G7Zu3Yrm5mYEBQXhscceg91ux9atW7Fq1SqMHTsWf/3rXxEZGenu8ok8AsOF6Aa5ubm499574e/vj61btyIuLg5CCJjNZkyePBk1NTUYMWIEMjMzERISAikljh07hoULFyIlJQXvvPMOfHx83P3PIHI7TugTXed0OrF27VpUV1dj06ZNarC0RwiB5ORkvPHGG/j888+xZ8+eHqqWyLMxXIiuy8/Px86dO/Hwww8jOTn5lsHSQgiBhx56CImJidiyZQscDoeLKyXyfFziQnRdeno66uvrMX/+fFy5cgUNDQ3qY4WFhXA6nQAAm82GCxcuICgoSH08KioKDz/8MFatWoWSkhJER0f3eP1EnoThQnTd5cuX4efnh5iYGDz11FM4fvy4+piUEs3NzQCA4uJi/OxnP1MfE0Lgj3/8I8aPHw+r1Yri4mKGC/V7DBei6xobG2EwGODj44Pm5mY0NTW1+Twp5U8eczgcMBqNrUKIqD9juFC/V1xcjNOnT+P8+fOwWq2wWCyYPn06/P391ec0NjYiPT1dDZGkpCR146QQAsOGDUNZWRkcDgfy8/MxdepU+Pr6uuufROR2XIpM/U5JSQkyMzNx+vRpZGZm4tq1axBCwN/fH8eOHcOmTZvwq1/9qtU1ZrMZU6dORU1NDYYPH47Tp08jODhYfVwIgRdffBFvvvkmFEWBr68vEhMTkZKSgpSUFEyfPp1LlKlfYbhQn1deXq4GyenTp1FUVAQAMJlMSEhIwJQpUzBp0iTY7XYkJycjJCQEe/bsaTVhf7N9LsAPw2TFxcWYNWsW7r//fjzxxBM4fPgwDh06hCNHjsBiscBoNOKOO+5Qw2bq1Knw9vZ2y9+DqCcwXKjPqaysbBUmV69eBQDExsaqYTJ58mQMGDDgJ9e+/fbb+Pd//3e8/PLLWLFihTr01V64NDU14fnnn8fOnTvxj3/8A6NGjVJ/n9PpxLlz53D48GEcPHgQR48eRW1tLfz8/JCUlITZs2cjJSUFCQkJ8PLy6oG/DlHPYLhQr1dVVYUzZ86oYXLlyhUAwIgRI5CQkICEhARMnjwZoaGht/xdDQ0NWLJkCXbv3o1XX30VS5cuha+vL7777jtMmzZNHRbLyMhAcHAw6urqsGbNGvz5z3/GW2+9hcWLF7f7+x0OB7KysnDo0CEcOnQIx44dQ319PQICAjBjxgw1bCZNmsTDMKlXY7hQr2OxWFqFidlsBgAMGzZM7ZkkJCQgLCysS7+/vLwczz77LL744gukpqZi+fLlGDNmDHJycqAoCry9vXH77bcjIyMDb775Js6ePYvVq1dj6dKl0Ov1nXotu92OzMxMNWyOHz8Oq9WKoKAgJCcnq2EzYcKETv9uIndiuJDHq62tVcMkMzMTeXl5AIDo6OhWYaLlqcQNDQ3YsmULNm7ciNLSUsTExMBkMiEwMBDV1dXIyclBcXExEhISsHLlSsyaNQs6XfcPvLDZbDh16pQ6jJaeno6mpiYEBwdj5syZatiMHz9ek9cjchWGC3mc+vp6ZGVlqfMmubm5kFJi8ODBapBMmTKlR04gLikpwYEDB3D48GGYzWY0NTUhJCQE48aNwz333IPp06fDz8/PZa/f3NyMr7/+Wg2bkydPorm5GaGhobjzzjvVsBk7dmyHj6sh6gkMF3I7q9WKrKwsdZirZfgpIiICU6ZMUf8bPHiwW+t0Op2QUkKn07mt19DU1ISTJ0+qw2gnT56E3W5HeHh4q7AZPXo0w4bciuFCPc5qteLcuXPqMNfFixehKArCw8PVCfgpU6ZgyJAh/IC8BavVihMnTqhhk5GRAYfDgcjISMyaNUsNG5PJxL8l9SiGC7lcU1MTsrOz1WGuCxcuwOl0IjQ0tNUw19ChQ/kB2E319fU4ceIEDh48iEOHDuH06dNwOp2IiopqFTYxMTH8W5NLMVxIczabDdnZ2WrP5Pz583A4HAgJCVF7JgkJCRg+fDg/4FystrYW6enpaticOXMGiqJg6NChrcJm+PDh7i6V+hiGC3VbyxH0LT2T8+fPw2azISgoqNUwF78tu5/FYsGxY8fUYbSzZ89CSonhw4e3CpuhQ4e6u1Tq5Rgu1Gl2ux0XL15Uw+TcuXOw2WwIDAzE5MmT1TCJjY3lclkPV1VVhWPHjqk9m+zsbABATEwMUlJS1LCJiopyc6XU2zBc6JYcDgcuXbqkruY6d+4cmpqa4O/vr4ZJQkICRo4cyTDp5SoqKnD06FE1bC5cuADgh3PYWsJm1qxZGDRokJsrJU/HcKGfcDqdyMnJwenTp9UwsVqt8PPzw8SJE9WeyahRo7hrvI8rKyvDkSNH1LC5fPkyAGD06NGtwkbLDazUNzBcCIqiIDc3Vx3mysrKQkNDA3x9fdUwSUhIwJgxY3jeVT937dq1VmHTclrC2LFj1bC58847u3z0DvUdDJd+SFEUfPvtt2rPJCsrC3V1dfD29saECRPU5cFxcXE8qZfaVVRUpJ4ecOjQIfWct/j4eDVsZs6cqZ4gTf0Hw6UfkFLCbDarYXLmzBnU1tbC29sb48ePV8Nk7NixvMcIdUtBQUGrsPn+++8hhMDEiRPVsElOTm7zdgfUtzBc+iApJa5cuaJOwJ85cwbV1dUwGAwYP368OswVHx/PMCGX+u6771qFTWFhIXQ6HSZPnqyGzYwZMxAYGOjuUkljDJc+QEqJq1evqj2TzMxMVFVVQa/XY+zYsWrPJD4+nvd1J7dp6UEfOnRIDZtr165Br9djypQpatgkJSXB39/f3eVSNzFceiEpJYqKitQgyczMRHl5OXQ6HeLi4lqFiStP7CXqDikl8vLy1LA5fPgwSktLYTAYMG3aNDVsEhMT+T7uhRguvcS1a9da9UxKS0uh0+kwatQo9dTgCRMm8Bsf9VpSSly+fLlV2FRUVMDb2xvTp09HSkoKUlJSkJiYyB54L8Bw6SUeeeQRFBQUYOTIkWrPZOLEiRyrpj5LURRcvHhRDZsjR46gqqoKn3zyCRYsWODu8ugWGC69RFNTE7y8vLhpkfqtGz+qeEad52O4EBGR5rjdWgMOhwMlJSVQFMXdpXSbEAKDBg3i5knqFLvdjoKCgj7TBoYNG8Zl+t3EcNFAaWkp1q9fjzFjxri7lG67dOkSVqxYgejoaHeXQr1IYWEhnnnmGSQkJAD4YQirtw5dZWZmYvPmzYiNjXV3Kb0aw0UjJpMJzzzzjMt+f0FBAfbt24eysjKMHTsWd911l0sm8zdt2gSOlFJnSSkRHx+P8ePH48svv0RiYiKWLl3q7rK6ZMWKFWwDGmC4aEzrb2uKouDIkSP4wx/+gPLycgDA3//+d+zevRurVq3C4MGDNXtNNijqrqNHj+Kjjz5CbW0tnnrqqV63AIVtQDu8+YYHk1IiOzsba9euRXl5uXofdKPRiLNnz2LdunVoaGhwd5lEqri4OACA2WxGU1OTm6shd2K4eLDa2lq89dZbsFgsuP3227Fx40asX78eK1asgK+vL06dOoVPP/2U37bIY7Tc46ekpARVVVXuLofciOGioTNnzuCjjz7CoUOHuv2BL6XEzp07cfnyZfj7++M//uM/MHToUOj1esyZMwcPPfQQFEXBtm3bcO3aNY3+BUTdM3z4cPj7+6O2thbff/+9u8vpECklKisrUV5eDofD4e5y+gyGi4a++uorpKWlYc+ePd0Ol8rKSmzfvh1SSsybNw/x8fHq3IrBYMDChQsRERGB8vJyfPHFF+y9kEeIiIhAREQE7HY7cnNze837cv369Zg0aRJWrFjRJ5ZTewKGi4Za7lFRU1PTrTeolBL79+/HtWvXEBISggULFvxkYjQyMhJz584FAOzZswc1NTVdL5xII/7+/hgxYgQA4Pz5826upmPsdjtOnTqFa9euwcvLq9cuofY0DBcNBQcHAwAaGhrgdDq7/HusVqvaG5k9e3abe06EEPj5z3+OwMBAFBcXIyMjo9d8S6S+q+U2DwBw8eLFbrWDnlJVVYW8vDwIITB16lSGi0YYLhpq6bk0NDTAbrd36Xe0rBAzm83w9fXFfffdd9M3+9ChQzFp0iQoioKvvvqqVzRk6vvi4+MBAN9++y3q6+vdXM2t5efno7KyEgEBARg3bpy7y+kzGC4aCgoKgk6ng9Vqhc1m69LvkFLiq6++gsPhwJgxY2AymW4aLnq9HnfffTeEEMjOzlb3wXT29err61FTU9PlmoluNHr0aHh7e6OsrAzFxcXuLqddUkpkZmbCbrcjOjoaQ4YMcXdJfQbDRUMt4dLc3IzGxsYu/Y7q6mpkZGQAAObMmdPu+UZCCEyePBkhISGwWCzIysrq0tDY22+/jYULF+Ivf/kLh9ao22677TaEhobCarXi8uXL7i6nXVJKnDx5EgB4cz2NMVw0FBAQAIPBAJvN1qXNjVJKnD17FhUVFQgKCsIdd9xxy/HfsLAwjBs3DlJKpKendzocpJQoKSlRl2FyvJm6KyQkBMOHD1ffz578haWhoQHZ2dkA0KH2Rh3HcNGQv78/vL294XA4UFdX1+nrpZQ4evQoFEVBXFwcBg0adMtr9Ho97rjjDgA/rM6pra3t1GsqiqLW2jJnRNQd3t7eGD9+PADg7NmzHr209/vvv8fVq1fh7e2NKVOmuLucPoXhoiFfX18YjUZIKbu0NLi2thZZWVkAgJkzZ3boXCYhBCZOnAg/Pz+Ul5fju+++69RrOhwOddK1ZbUbUXe1fFBfvny5S1+0eoKUEllZWbBarYiMjITJZHJ3SX0Kw0VDPj4+8Pf3h5QS1dXVnb4+Ly8PZWVl8PPzQ0JCQoe76FFRURgyZAgcDken511sNhusViuEEOy5kCaEEIiPj4ePjw+uXbuGq1evurukmzp69CiklBg/fjxCQkLcXU6fwnDRkMFgQFBQEIAf1s535kNeSomvv/4aTqcTw4cPR1RUVIev9fX1VYchsrKyOrUkuampCY2NjRBCqLUTddeIESMQHh6OxsZGnD9/3iPnXaxWK06dOgUAmDFjBnQ6fhxqiX9NDel0OvXbf2cP7bPZbDhz5gwAICEhAT4+Pp26ftKkSRBCwGw2d2pIrmXZ9I3BSNRdISEhGDNmTKvVWJ6moKAAZrMZXl5eSEpK4mS+xhguGhJCIDQ0FMAPS4o7822trKwMV65cgU6nw5QpUzr1RhdCYNSoUTAajaiurkZBQUGHr62rq4PdboeXlxcCAgI6fB1Re/R6PaZOnQoAOH36tMftoZJS4tSpU2hoaEBkZGSfuIusp2G4aOzGcOnoKhkpJS5evIiGhgaEhoZi5MiRnX7dyMhIDB48GA6HAxcuXOhwsNXW1kJRFHUxApEWhBBITEyETqdDfn4+SkpK3F1SK1JK/OMf/4CUEpMmTVLbLWmH4aKxgQMHAvjh8MrOHN995swZSClx++23d2li3cfHB6NGjQIAfPPNNx0Ol5YQ9PPz6/RQHFF7WibJq6ur1b0knqKurg5ff/01ACAlJYXzLS7Av6jGQkNDIYRAfX19h+/E19zcjG+++QYAMHHixC7dGlYIoZ6LlJ+fD6vVestrpJTq3FBQUBC8vLw6/bpENxMZGYnRo0dDURR1VZanyMnJwffffw9fX1/MnDmT8y0uwHDRWEhIiHq+WEc+4AGgtLQUxcXF0Ov1mDBhQpfe6C3zLl5eXigvL0dZWVmHrmsJlwEDBvS6+52TZ/P29saMGTMAAMePH/eYeRcpJQ4dOoTm5maMGDGiS8PQdGsMF40FBwfD29sbNputQ7vlpZTIycmB1WpVj83oqujoaAQHB6OpqQn5+fkd+qZYWVkJ4P+HIpFWhBCYNWsW9Ho9Ll++7DH7XRwOB/bv3w8ASE5O5kIWF+GnicYCAgLg6+sLh8PR4Y2U586dg5QSMTEx3drIGBgYiGHDhkFKiQsXLtzy+YqiqDWGhYV1+XWJbmbixImIiIhATU0NTp486RFDY4WFhcjKyoJOp0Nqaqq7y+mzGC4aMxqNCAgIgKIoqKiouOXz7XY7Ll68CAAYN25ct4am9Ho94uLiAPwwpnyrzZQ3BmB4eHiXX5foZsLDwzF16lRIKfHll1+6PVyklDh27Biqq6sRERGB6dOnc77FRRguGvPx8VGPkSgtLb1lY6qsrMTVq1eh0+kQHx/frTe6EAJjxoyBEAIFBQW33EzZ3NyMuro6CCHYcyGXuLF3kJ6erg7DuouiKNi5cyeklEhMTERkZKRb6+nLGC4a0+l06nLk0tLSWz4/Pz8f9fX1CAwMRExMTLdfPzY2FkajETU1NSgqKmr3uS2LDnQ6HcOFXEIIgZSUFAQFBaGoqMjtt+MuKSlBeno6hBC4//77Oc/oQvzLakwIoX4bKi8vb3cjpZQS58+fh6IoGDp0qCYbuSIiIjBw4EDY7Xbk5ua225Bra2vR1NQELy8vnohMLhMTE4NJkybB6XTis88+c1u4SClx5MgRlJaWIiwsDCkpKRwScyGGiwvcGC52u/2mz3M6nTh//jyAH+ZbtNhnYjQaERsbCwC3nNSvrKyEw+GA0WhEYGBgt1+bqC1eXl546KGHAAD79+/v0FykKzidTuzYsQOKoiApKQnR0dFuqaO/YLhoTAiBwYMHQwgBi8XS7kbKmpoaXLlyRT2iXItvUUIIdVI/Ly8Pzc3NN31uWVkZFEVBYGAgb+9KLiOEwM9//nOEhISgsLBQPXalp129ehVHjx6FEALz58/nvi4XY7i4QEREBAwGAxoaGtqdVP/uu+9gsVjg5+en2Y2KhBAYPXo0dDodSkpKbno6s5RSnRMKDQ2Ft7e3Jq9P1JYRI0Zg1qxZUBQFH374Ybs9eleQUmLXrl2oqKhAVFQU7rrrLg6JuRjDxQVCQ0NhNBrR1NSE8vLyNp/TMt/idDoRFRWl6aqV2267DYGBgaivr8eVK1du+rxr164B+CEM+S2OXEmv1+OXv/wlDAYDjh492uP3eGlsbMQnn3wCAEhNTe3QLcSpexguLhAYGIgBAwZAURQUFxe3+Ryn06ke5hcXF6fpoZEhISGIioqCoig3PSHZ6XSqPZfBgwdr9tpEbRFCYPbs2YiLi0N9fT22bNnS4VPDu6vleP2srCz4+Phg4cKF7LX0AIaLC/j6+iIiIgLAD+O8bX2419XVIS8vD0IITJw4UdPX9/b2xujRowEAFy9ebLMR22w2dWJ1yJAhmr4+UVuCgoLw5JNPQgiBHTt2ICcnp0d6L4qi4N1330VzczMmT57MjZM9hOHiAnq9Xl2JUlBQ0GYD+u6771BVVQWj0Yi4uDjN3+wttz02m82or6//yeP19fWorq6GTqdTFyAQuZIQAgsWLIDJZEJVVRXeeuutTt2Su6tyc3Oxe/duCCGwaNEi3reohzBcXKTlAMrCwsKfnAYrpcTZs2fhcDgQFRWl+fhvywnJvr6+qKysRGFh4U+eU1FRgcbGRvj4+Ki9LCJXCw8Px7Jly6DT6bBt2zaXH8WvKAreeecdVFVVITY2Fg8++CC/SPUQhosLCCEwYsQI6HQ6lJeX/+R0ZKfTiczMTABAfHw8fH19Na9h8ODBCA8Ph81ma3Pepbi4GHa7HYGBgepxNUSuJoTA448/jsTERDQ0NODFF19EVVWVSwJGSolvv/0WH3/8MYQQWLJkiXp6Brkew8VFhgwZAl9fX9TX16urslpUVFQgLy8POp1Ovc+41vz8/NR5l7Nnz7ZqvFJKXLlyBVJKREREwN/f3yU1ELUlMDAQr732GoKCgnDq1CmsWrXKJUuTFUXBxo0bUVZWhuHDh+OXv/wley09iOHiIgMHDsTAgQPhcDha3VtFSolvvvkGtbW1CAoKwtixY13yhhdCYNKkSQCAy5cvo6GhQX1MSgmz2QwAGDZsGAwGg+avT3QzQggkJyfjN7/5DXQ6HbZs2YJNmzZ16rbgt9KyQuzDDz+EEAK//vWvuSqyhzFcXMRoNKoHUbbcwhj44U1//PhxKIqCUaNGuayb3rLr38fHB2VlZa32u9hsNvV/a7V5k6gz9Ho9nn/+eSxYsAAOhwOvvPIK/vSnP2l2t8qGhga88sorqK2txaRJk/DEE0+w19LDGC4uIoRQV2xdunQJjY2NAACLxaLOtyQnJ7t082J0dDSGDBkCm82G06dPq72nqqoqlJWVQafTwWQysdGRW/j5+SEtLQ3z5s1DU1MTXnzxRfz2t7+FxWLp1hyMoih4++23cfDgQRiNRqxatYoHs7oBw8VFWnoOXl5eKC4uRmFhIaSUOHPmDMrKyhAYGIjExESXfrAbjUZMnjwZAHDy5EnY7XZ1SKzlmP9hw4a57PWJ2iOEQGhoKN599121B5OWloZHHnmkyzv4pZTYs2cP1q9fDykllixZgtTUVH6BcgOGiwvFxMRg0KBBaGpqQnp6Oux2O7744gsoioKJEyf2yObFGTNmQK/XIycnR12SnJWVBUVRMGzYME2O+SfqqpaA2bJlC1566SX4+fnh4MGDuPfee7F582Y0NDR0OGSklEhPT8fTTz+N2tpa3HHHHXjllVc4p+gmDBcXCggIQFJSEgBg7969OHjwIE6fPg29Xo8HHnjA5ed5CSEwbtw4DB48GFarFQcPHkRzczNOnToFAJg0aZImx/wTdYcQAv7+/vjd736HTz75BKNHj0ZpaSmWL1+OBx98EIcPH1Z73TfjdDqxc+dO/PM//zOKiopgMpnw5z//mTfBcyOGiwsJITBv3jz4+/vDbDZj7dq1sNlsmDhxIqZNm9YjXfWgoCCkpKQAAHbv3o1jx47BbDbDy8sLM2bM4HABeQy9Xo+5c+di3759+NWvfgVvb28cPHgQ9913Hx599FHs2rULVVVVcDqdkFJCSgm73Y5Lly7h3/7t3/CLX/xCDZYPP/xQveU3uQf7iy5mMpnw2GOP4b333kNjYyMiIyPx3HPPuWTjZFtabue6c+dOFBYWYv369bDZbBg3bpy6D4bIUwghEBUVhbfffhuPPPII1qxZgxMnTuDzzz/Hrl27EB0djbFjx2LIkCFwOBzIy8tDdnY2ampqoNPpcNddd+FPf/oTRo0axWBxM4aLi+n1eixevBgmkwlFRUVITk7G8OHDe/SNf9ttt2HhwoX4y1/+gtraWhiNRjzxxBM9FnBEnSGEgMFgwJw5c5CUlIR9+/bhnXfeQXp6Oq5cufKT20jo9XqMHDkSS5cuxaJFixAYGMhg8QAMF421NS7s5eWF2bNn3/J5riKEwMKFCxEaGors7GzMnDkTycnJPV4H9Q9avqf8/Pzw0EMP4b777oPZbMbXX3+Nc+fOobS0FHq9HrfddhsSExORmJjY6hgjvq/dj+Gikfz8fGzZssXdZdxSZGQkcnNzkZub2+bjubm5/NZHnSaEwPnz5/Hqq6+6/LUCAwMREBCgvk8zMjKQkZGh2e8/d+4c24AGhGTEd5vdbkdRUVGP3fzIlYQQiI6O5ioy6hSbzQaz2dwjR+i7mk6nQ2xsLG/93U0MFyIi0hyXIvcSiqKgsbGxT/SOiLqqZQkyvxN7PoZLL5Gbm4s777zzpnMlRP1BVlYWDAYDsrKy3F0K3QLDhYiINMdwISIizTFciIhIcwwXIiLSHMOFiIg0x3AhIiLNMVyIiEhzDBciItIcw4WIiDTHcCEiIs0xXIiISHMMFyIi0hzDhYiINMdwISIizTFciIhIcwyXXkBKierqathsNlRXV/NGSdQvtbQDAGwHvQDDxYNZLBakpaXBZDIhKSkJ2dnZSEpKgslkQlpaGiwWi7tLJHK5G9vBnDlzoCgK5syZw3bg4YRk/HukvXv3Yv78+bBarQDQ6luaEAIA4Ofnhx07diA1NdUtNRK5GttB78Vw8UB79+7FvHnzIKWEoig3fZ5Op4MQArt27WLDoj6H7aB3Y7h4GIvFgujoaDQ2NrbboFrodDoYjUYUFhYiODjY9QUS9QC2g96Pcy4e5v3334fVau1QgwIARVFgtVrxwQcfuLgyop7DdtD7sefiQaSUMJlMMJvNnVoJI4RATEwM8vLy1HFoot6K7aBvYLh4kIqKCoSHh3fr+rCwMA0rIup5bAd9A4fFPEh9fX23rq+rq9OoEiL3YTvoGxguHiQgIKBb1wcGBmpUCZH7sB30DQwXDxIWFobY2NhOjxcLIRAbG4vQ0FAXVUbUc9gO+gaGiwcRQuC5557r0rXLli3jJCb1CWwHfQMn9D0M1/cTsR30Bey5eJjg4GDs2LEDQgjodO3/39OyM/nTTz9lg6I+he2g92O4eKDU1FTs2rULRqMRQoifdPNbfmY0GrF7927cc889bqqUyHXYDno3houHSk1NRWFhITZs2ICYmJhWj8XExGDDhg0oKipig6I+je2g9+KcSy8gpURVVRXq6uoQGBiI0NBQTlpSv8N20LswXIiISHMcFiMiIs0xXIiISHMMFyIi0hzDhYiINMdwISIizTFciIhIcwwXIiLSHMOFiIg0x3AhIiLNMVyIiEhzDBciItIcw4WIiDTHcCEiIs0xXIiISHP/D4srzhLUwrHqAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 500x400 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "187d19f9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "saving model version 0.2\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "| train_loss: 1.54e-02 | test_loss: 4.73e-02 | reg: 7.50e+00 | : 100%|█| 20/20 [00:02<00:00,  6.93it"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "saving model version 0.3\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "model = model.refine(20)\n",
    "model.fit(dataset, opt=\"LBFGS\", steps=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8d50bcef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxPElEQVR4nO3deXSTZd4+8OtO0iVtWtKNpS1LUyqLgkgpSym0LEMd0VHBceHMjLjr4IZzjjrMT0XG5dWRERAdfdHXQXQOKOVVthFGhAJlKVRAdihla0v3plvaNM1z//4oed4WAUGe9Enb63MOR2mS5tuSO1fu9RFSSgkiIiINGfQugIiIOh6GCxERaY7hQkREmmO4EBGR5hguRESkOYYLERFpjuFCRESaY7gQEZHmGC5ERKQ5hgsREWmO4UJERJpjuBARkeYYLkREpDmGCxERaY7hQkREmjPpXQBReyClRHl5OWpra2GxWBAREQEhhN5lEfks9lyILsNut2P+/PlISEhAVFQU4uLiEBUVhYSEBMyfPx92u13vEol8kuCVKIkubt26dZg6dSocDgeA5t6Lh6fXEhQUhIyMDKSnp+tSI5GvYrgQXcS6deswefJkSCmhKMol72cwGCCEwJo1axgwRC0wXIguYLfbERsbi/r6+ssGi4fBYIDZbEZ+fj6sVqv3CyRqBzjnQnSBxYsXw+FwXFGwAICiKHA4HPjss8+8XBlR+8GeC1ELUkokJCQgLy8PV9M0hBCw2Ww4fvw4V5ERgeFC1EpZWRmioqKu6fEREREaVkTUPnFYjKiF2traa3p8TU2NRpUQtW8MF6IWLBbLNT0+JCREo0qI2jeGC1ELERERiI+Pv+p5EyEE4uPjER4e7qXKiNoXhgtRC0IIPPXUU7/osU8//TQn84nO44Q+0QW4z4Xo2rHnQnQBq9WKjIwMCCFgMFy+iXh26K9YsYLBQtQCw4XoItLT07FmzRqYzWYIIX4y3OX5mtlsxtq1azFp0iSdKiXyTQwXoktIT09Hfn4+5s2bB5vN1uo2m82GefPmoaCggMFCdBGccyG6AlJKbNy4ERMmTMCGDRswbtw4Tt4TXQZ7LkRXQAihzqlYrVYGC9HPYLgQEZHmGC5ERKQ5hgsREWmO4UJERJpjuBARkeYYLkREpDmGCxERaY7hQkREmmO4EBGR5hguRESkOYYLERFpjuFCRESaY7gQEZHmGC5ERKQ5hgsREWmO4UJERJpjuBD9DJfLhYKCAhw+fBgAcOLECVRUVEBRFJ0rI/JdvMwx0SXY7XZkZGTgiy++wMGDB1FTU4PGxkYEBgYiKioKY8aMwUMPPYTRo0fDZDLpXS6RT2G4EF3E9u3bMXPmTPz4449ISkrC5MmTMXjwYFgsFtjtduTk5GDVqlXIzc3FPffcg9deew1RUVF6l03kMxguRBdYv349pk+fDovFgjfffBO33HILGhsbsXTpUjidToSGhuLee++Fy+XC0qVLMXv2bFx//fVYsmQJunXrpnf5RD6B4ULUwrFjx3DzzTcjODgYS5cuxcCBAyGEQF5eHoYOHYqqqirExcUhJycHYWFhkFJi69atmDZtGtLS0vDxxx8jICBA7x+DSHec0Cc6z+1244033kBlZSUWLlyoBsvlCCGQkpKCt99+G9988w2+/fbbNqqWyLcxXIjOy83NxapVqzBlyhSkpKT8bLB4CCFwxx13YOTIkVi0aBGampq8XCmR7+MSF6Lztm3bhtraWkydOhWnTp1CXV2delt+fj7cbjcAoLGxEQcPHkRoaKh6e3R0NKZMmYLZs2ejqKgIsbGxbV4/kS9huBCdd+TIEQQFBcFms+Gxxx5DVlaWepuUEk6nEwBQWFiIX/3qV+ptQgjMnTsXgwYNgsPhQGFhIcOFOj2GC9F59fX1MJlMCAgIgNPpRENDw0XvJ6X8yW1NTU0wm82tQoioM2O4EJ3XtWtX1NfXw263Y8SIEQgODlZvq6+vx7Zt29QQSU5OVjdOCiHQq1cvlJSUwGAwICwsTK8fgchnMFyIzktMTITL5UJ2djbeeuutVrfl5eUhKSkJVVVV6NatG5YtWwar1areLoTArFmz0L17dw6JEYGrxYhUw4cPh81mw+LFi1FXVwej0djqj4cQAgaDQf26wWDAuXPnsHz5ckyePBldunTR8acg8g0MF6LzIiIi8OSTT+KHH37AggULrnhJsdPpxF//+lfU19fjscceu+IlzEQdGYfFiFqYPn06Nm/ejLfeegtBQUF44oknEBgYCAAwmUwwmUxqL0ZKiZqaGrz++utYunQp3n33XfTr10/P8ol8Bo9/IbpAaWkpZsyYgdWrVyM9PR0zZ87EgAEDcPToUSiKAn9/f/Tt2xfZ2dl45513sHfvXsyZMwdPPPFEq+Ezos6M4UJ0EXV1dVi0aBEWLFiA4uJi2Gw2JCQkICQkBJWVlTh69CgKCwuRmJiIV155BampqTAYOMpM5MFwIbqMoqIibNiwAZmZmdi3bx+ys7MxZswYjB49GpMmTcKIESMQFBSkd5lEPofhQnSFdu3aheHDh2PXrl0YNmyY3uUQ+TT244mukNFoVJchE9HlsZUQEZHmGC5ERKQ5hgsREWmO4UJERJpjuBARkeYYLkREpDmGCxERaY7hQkREmmO4EBGR5hguRESkOYYLERFpjuFCRESaY7gQEZHmGC5ERKQ5Xs+F6ApJKaEoCgwGA4QQepdD5NPYcyG6CryWC9GVMeldAJEWXC4Xzpw5A0VR9C7lmgkh0KtXL/j7++tdCtEvxnChDiE/Px9//OMfkZiYqHcp1ywnJwcffPAB4uPj9S6F6BdjuFCHIKXE4MGD8frrr+tdylU7cuQI/vu//xt33XUXkpOT8eKLL4JTodTeMVyow2lPk+35+fm47777sH//fuzcuRP/+c9/9C6JSBOcnSTSiZQSH374Ifbv3w8A2LNnD7Kzs3WuikgbDBcinZw+fRqLFy8G0LwKrbGxEVu3btW5KiJtMFyIdCClxCeffILCwkLExMRg3LhxMBqNyMvL43wLdQiccyHSwZkzZ/DPf/4TAPDAAw/gzjvvRFVVFRISEjBv3jxdayPSAsOFqI0pioJFixapvZaHHnoIvXr1AtDco2lPCxKILoXDYkRt7PTp02qv5cEHH0TPnj31LYjICxguRG1IURR89NFHOHfuHHr27ImHHnqIPRXqkBguRG0oLy8Pn332GQDg4YcfRmxsrM4VEXkHw4WojSiKgg8//BDFxcXo3bs3HnjgAfZaqMNiuBC1ASkljh07hs8//xxCCDz66KOIjo7Wuywir2G4ELUBRVGwYMEClJaWwmaz4f7772evhTo0hguRl0kpsX//fixbtgxCCMyYMQPdu3fXuywir2K4EHlZU1MT5s6dC7vdjoEDB+J3v/sdey3U4TFciLxISont27fj66+/htFoxHPPPYfw8HC9yyLyOoYLkRc1NDTgv/7rv+BwODB8+HBMnTqVvRbqFBguRF4ipcTq1avx/fffIyAgAC+88AIsFoveZRG1CYYLkZdUVFTgjTfegMvlwq9//WtMmjSJvRbqNBguRF7gOeZl//79sFqtmDVrFvz9/fUui6jNMFyINCalxOHDh/Hee+9BSokHH3wQQ4YMYa+FOhWGC5HGnE4nXn31VZSUlCAhIQEzZ86E0WjUuyyiNsVwIdKQlBIZGRlYuXIl/Pz88PLLL6NHjx56l0XU5hguRBqRUuLUqVN45ZVX4HK58Jvf/AZTpkzhcBh1SgwXIo04nU7MmjULJ0+eRExMDObMmYOAgAC9yyLSBcOFSAOKouCTTz7BihUr4Ofnh1dffRX9+vVjr4U6LYYL0TWSUmLz5s2YPXs2mpqacO+992LatGkMFurUGC5E10BKiePHj+OJJ55ARUUFbrrpJrz55pvc00KdHsOF6BeSUqK4uBgPP/wwjh07hh49euAf//gHunfvzl4LdXoMF6JfQEoJu92Oxx9/HFlZWQgNDcV7772HYcOGMViIwHAhumpSStTW1uKpp57C6tWrERgYiLfffhu33347g4XoPIYL0VWQUqKurg4zZ87EsmXL4Ofnh1deeQUPPPAADAY2JyIPtgaiKySlRE1NDZ555hksXrwYBoMBzz//PJ599lmYTCa9yyPyKWwRRFdASomysjI8+eSTyMjIgMlkwp/+9Ceedkx0CQwXop8hpcSJEyfw2GOPYdOmTQgICMCLL76IF154gTvwiS6B4UJ0GYqi4LvvvsPTTz+N48ePw2Kx4LXXXsPjjz8OPz8/vcsj8lkMF6KLkFKioaEB77//Pl5//XVUV1cjOjoa8+fPx+23384j9Il+BsOF6AJSSuTl5eEvf/kLVqxYAUVRkJycjIULF2Lw4MFcbkx0BbhajOg8KSWcTif+53/+B+PHj8dXX30Fo9GIRx55BP/7v//LYCG6Cuy5UKcnpYTb7cbu3bvx97//HStXroTL5ULv3r3xyiuv4L777oOfnx+DhegqMFyoU1MUBSdOnMDcuXPxr3/9C3V1dfDz88O0adPw17/+Fb1792aoEP0CDBfqdKSU6lUjP/roI3z22WcoKSmBwWDA8OHD8eyzz+KOO+6Av78/g4XoF2K4UKchpYTL5cLevXvx6aef4uuvv0ZJSQmEEIiPj8fMmTPxu9/9DhaLhaFCdI0YLtSheeZT8vPz8d1332HZsmXYsWMHHA4HhBDo06cPHn30UUyfPh1du3ZlqBBphOFCHY6UEoqi4Ny5c/juu++QkZGBnTt3oqKiAlJK+Pn5ITExEQ888ADuvPNOdOvWjaFCpDGGC3Uoy5cvx+bNm3Hq1Cns3r0bpaWlkFLCYDAgOjoaEydOxL333ovk5GQEBwczVIi8hOFCHcrmzZvxwQcfAACEEOjevTvGjRuHKVOmIDk5GVFRURBCMFSIvIzhQh1KYmIifvvb36Jbt24YPXo0UlJSLnrZYSmlThUSdQ4MF+oQhBDYv38/goKC0K9fPwghcOjQIRw6dEjv0q7avn372LOidk9IfoSjDqCxsRF5eXlwu916l3LNDAYD4uPjeZ0YatcYLkREpDkOixFdoZafwzhsRXR5PBWZ6Art2bMHRqMRe/bs0bsUIp/HcCEiIs0xXIiISHMMFyIi0hzDhYiINMdwISIizTFciIhIcwwXIiLSHMOFiIg0x3AhIiLNMVyIiEhzDBciItIcw4WIiDTHcCEiIs0xXIiISHMMFyIi0hzDhegKSClRWVkJAKisrAQv4Ep0eQwXosuw2+2YP38+EhISMHHiREgpMXHiRCQkJGD+/Pmw2+16l0jkk4TkRzCii1q3bh2mTp0Kh8MB4OKXOQ4KCkJGRgbS09N1qZHIVzFciC5i3bp1mDx5MqSUUBTlkvczGAwQQmDNmjUMGKIWGC5EF7Db7YiNjUV9ff1lg8XDYDDAbDYjPz8fVqvV+wUStQOccyG6wOLFi+FwOK4oWABAURQ4HA589tlnXq6MqP1gz4WoBSklEhISkJeXd1UrwoQQsNlsOH78uDofQ9SZMVyIWigrK0NUVNQ1PT4iIkLDiojaJw6LEbVQW1t7TY+vqanRqBKi9o3hQtSCxWK5pseHhIRoVAlR+8ZwIWohIiIC8fHxVz1vIoRAfHw8wsPDvVQZUfvCcCFqQQiBp5566hc99umnn+ZkPtF5nNAnugD3uRBdO/ZciC5gtVqRkZEBIQQMhss3Ec8O/RUrVjBYiFpguBBdRHp6OtasWQOz2QwhxE+GuzxfM5vNWLt2LSZNmqRTpUS+ieFCdAnp6enIz8/HvHnzYLPZWt1ms9kwb948FBQUMFiILoJzLkRXQEqJjRs3YsKECdiwYQPGjRvHyXuiy2DPhegKCCHUORWr1cpgIfoZDBciItIcw4WIiDTHcCEiIs0xXIiISHMMFyIi0hzDhYiINMdwISIizTFciIhIcwwXIiLSHMOFiIg0x3AhIiLNMVyIiEhzDBciItIcw4WIiDTHcCEiIs0xXIiISHMMF6Kf4XK5UFBQgMOHDwMATpw4gYqKCiiKonNlRL6LlzkmugS73Y6MjAx88cUXOHjwIGpqatDY2IjAwEBERUVhzJgxeOihhzB69GiYTCa9yyXyKQwXoovYvn07Zs6ciR9//BFJSUmYPHkyBg8eDIvFArvdjpycHKxatQq5ubm455578NprryEqKkrvsol8BsOF6ALr16/H9OnTYbFY8Oabb+KWW25BY2Mjli5dCqfTidDQUNx7771wuVxYunQpZs+ejeuvvx5LlixBt27d9C6fyCcwXIhaOHbsGG6++WYEBwdj6dKlGDhwIIQQyMvLw9ChQ1FVVYW4uDjk5OQgLCwMUkps3boV06ZNQ1paGj7++GMEBATo/WMQ6Y4T+kTnud1uvPHGG6isrMTChQvVYLkcIQRSUlLw9ttv45tvvsG3337bRtUS+TaGC9F5ubm5WLVqFaZMmYKUlJSfDRYPIQTuuOMOjBw5EosWLUJTU5OXKyXyfVziQnTetm3bUFtbi6lTp+LUqVOoq6tTb8vPz4fb7QYANDY24uDBgwgNDVVvj46OxpQpUzB79mwUFRUhNja2zesn8iUMF6Lzjhw5gqCgINhsNjz22GPIyspSb5NSwul0AgAKCwvxq1/9Sr1NCIG5c+di0KBBcDgcKCwsZLhQp8dwITqvvr4eJpMJAQEBcDqdaGhouOj9pJQ/ua2pqQlms7lVCBF1ZgwX6vROnTqFTZs2YevWrXA4HLDb7RgxYgSCg4PV+9TX12Pbtm1qiCQnJ6sbJ4UQ6NWrF0pKStDU1ITc3FwkJSUhMDBQrx+JSHdcikydztmzZ5GZmYnMzExs2rQJZ86cgRACvXv3xokTJ/D+++/j4YcfbvWYvLw8JCUloaqqCn369MHu3bthtVrV24UQmDVrFt555x0YjUYEBgZixIgRGDt2LNLS0pCUlMQlytSpMFyowyssLMSmTZuwefNmbNq0CSdPngQADB48WH3zT0lJgaIoSElJQVhYGL799ttWE/aX2ucCNA+TFRYWIjU1FbfddhumT5+OzZs3Y/PmzdiyZQvsdjvMZjNGjhypPl9iYiL8/f11+X0QtQWGC3U4xcXFrcIkNzcXAHD99derb+5jxoxBRETETx77/vvv409/+hP+3//7f3jxxRfVoa/LhUtDQwOeffZZrFq1Ct9//z369eunfj+3240ff/wRmzdvRmZmJrZu3Yrq6moEBQVh1KhRSE1NRVpaGm666Sb4+fm1wW+HqG0wXKjdKy0tRWZmphomR48eBQD0798fqampSE1NxdixY6/o7K+6ujo8+OCDWLt2LV599VU88cQTCAwMxMmTJzF8+HB1WCw7OxtWqxU1NTV4/fXX8dFHH+Hdd9/FAw88cNnv39TUhL1796r1ZmVloba2FhaLBcnJyWq9Q4YM4WGY1K4xXKjdKSsrw5YtW9R5k0OHDgEAEhIS1Dfn1NTUX3zOV2lpKWbMmIHVq1cjPT0dM2fOxIABA3D06FEoigJ/f3/07dsX2dnZeOedd7B3717MmTMHTzzxBIxG41U9l8vlwp49e9Se1rZt2+BwOBAaGorRo0erP8vgwYOv+nsT6YnhQj6vsrISW7ZsUd+A9+/fDwCw2WytwiQ6Olqz56yrq8OiRYuwYMECFBcXw2azISEhASEhIaisrMTRo0dRWFiIxMREvPLKK0hNTYXBcO0HXjQ2NiInJ0cNzu3bt6OhoQFWqxUpKSnqz3rDDTdo8nxE3sJwIZ9TVVWFrVu3qmGyb98+SCnRu3dvpKWlqfMmbbFRsaioCBs2bEBmZiby8vLQ0NCAsLAw3HDDDZg0aRJGjBiBoKAgrz2/0+nErl271LDZuXMnnE4nwsPDMWbMGDVsruQcNKK2xHAh3dXU1CArK0t9A92zZw8URUFMTAzS0tLUSe/evXvrWqfb7YaUEgaDQbdeQ0NDA3bu3KnOL2VnZ8PlciEyMhJjx45Vw6Zfv34MG9IVw4XaXG1tLbZv3672THJycuB2u9GjR49WPZO4uDi+Qf4Mh8OBHTt2qGGze/duNDU1oWvXrq2GDPv27cvfJbUphgt5necN0BMmu3btUt8AW4YJ3wCvXW1tLXbs2KFuEP3hhx/U4G4ZNgxu8jaGC2nOM3TjCZOdO3fC5XIhKioKY8eOVcOEQzfeV11dje3bt/9kyDE2Nlb9d0hNTdV9yJE6HoYLXTOn04ns7OyLTjq3nAfgpLP+7HY7srKy1E2dLRdLtAwbnupM14rhQletsbERu3fvVodeduzYoS6XbbmCictlfV9FRUWrxRSeZd5xcXHqv2NaWhp69Oihc6XU3jBc6Ge5XC788MMPaphs375d3eiXkpKizptwo1/7V15eftENqn379tVkgyp1HgwX+omWR5Rs2rRJvUKjxWJptWucR5R0fKWlpeoQWmZmpnq0Tr9+/VqFTWRkpM6Vkq9huJB6uKInTLKystTDFZOTk9WeydChQ3m4YidXVFTUKmw8h4IOHDiw1Tlu4eHhOldKemO4dEKKouDAgQNqmGzduhV2ux2BgYEYNWqUGibDhg3jsfB0WYWFhWrQZGZmqpczGDRokBo2nssYUOfCcOkEpJQ4dOiQGiZbtmxBRUUFAgICMGLECDVMhg8fzgta0TU5e/asuqEzMzNTvRDbjTfeqIbN6NGj0aVLF71LJS9juHRAUkocPXpUbeBbtmxBaWkp/Pz8MGLECLWRjxgxgpfiJa86depUq7ApKCiAwWDATTfdpL4Ok5OTERISoneppDGGSwcgpURubq7agDMzM1FSUgKTyYSkpCS1EY8cOdKrhywSXY6UEidPnmx1iemioiIYjUYMHTpU3WMzatQoBAcH610uXSOGSzvkaaQtw+TcuXMwGo1ITExU9yaMHDkSFotF73KJLsrzoajlnI3nQ9GwYcPUsPH2ydPkHQyXduL06dOtwiQ/P7/V8EJaWhpGjRrV6rrvRO2JZzi3ZdiUl5fD398fSUlJ6twgh3PbB4ZLOzFo0CAcP35cnRhNS0tDcnIyrFar3qUReYWiKDh8+LAaNJ6FKEuWLMFvf/tbvcujn8FwaScURYEQgmdzUaclpYSUku2gnWC4EBGR5nh2hwZcLhfOnj0LRVH0LuWaCSHQs2dPbp6kq8I2QBdiuGigoKAATz75JBITE/Uu5Zrl5ORg4cKFsNlsepdC7UjLNuAZvmqvJ2KzDWiD4aIBKSUGDx6MOXPmeO05jh49io8//hi33XYbxo4d67Xn+ctf/gKOlNLVklJi0KBBGDJkCFauXImxY8fiwQcf1LusK7Z+/XocPHgQ48ePR1NTE9uABhguGvPGRKPb7cbs2bPx9ddfY+PGjdi4caNXlhyzQdG1EELgu+++w7Jly2C323H//fe3i1OzFUXBkiVLsHz5ctx9992Ijo7Wu6QOoX32WzuZqqoq/PDDDwCaj9MoKCjQuSKiixszZgyEEDhw4ADKy8v1LueK2O127Nq1C0Bz/e11OM/X8LfYDpSUlKgNtb6+HmfPntW5IqKLS0xMhMViQUlJCY4cOaJ3OVfk0KFDOHfuHIKDgzFq1Ci9y+kwGC7tQFFRERoaGgA0d+HZcyFf1bNnT9hsNrhcLmRlZfn8UKuUEpmZmWhsbER8fDzi4uL0LqnDYLi0AwUFBXC73erfCwsLfb7RUudkNpuRlJQEAMjKymr1uvVFLpcLmzZtAgCMHj2aZ5hpiOHi46SUKCwsbPW14uJinaoh+nljx46FEAIHDx5EaWmp3uVcVmFhIQ4dOgSDwYDx48frXU6HwnBpB4qKilr9vaysjD0X8lmJiYno0qULSktLsX//fr3LuSQpJXbt2oWKigpERkbipptu0rukDoXh0g5c+OmvvLy8Q+yEpo4pJiYG1113HdxuNzIzM332g5CUEt999x2klLjxxhvRvXt3vUvqUBguPs7tdqOsrAwA1GuzVFVVoampSc+yiC4pICAAKSkpAIAtW7bA6XTqXNHF1dTUYPv27QCACRMmwGg06lxRx8Jw8XFNTU2w2+0AgL59+wIAamtr4XK5dKyK6NKEEBg3bhxMJhOOHDmCM2fO6F3SRR05cgSnT5+G2WxW54lIOwwXH9fY2IiqqioAQHx8PIQQcDgcPvtpkAiAOsxUU1ODbdu2+dzQmJQSGzduhNPphM1mw3XXXad3SR0Ow8XHNTQ0oK6uDgaDATabDUajEfX19aivr9e7NKJLioiIwPDhwyGlxLp163xujtDlcmHDhg0AmnflBwcH61xRx8Nw8XEOhwMNDQ0wGAzo2bMnTCYTnE4n6urq9C6N6JIMBgNuvvlmCCGwc+dOn1uSXFBQgAMHDsBoNGLixIl6l9MhMVx8XG1tLZxOJ/z8/NCjRw/4+fmhqakJtbW1epdGdElCCKSkpMBqtaKoqAjZ2dk+MzQmpURWVhbsdju6du2KxMREzrd4AcPFx1VXV8PlciEgIABRUVEICAiA2+1GTU2N3qURXVbPnj0xZMgQKIqCtWvX+lS4rF+/HlJKDBs2DF27dtW7pA6J4eLjqquroSgKAgMDERYWBrPZDEVRUF1drXdpRJdlMpnw61//GgCQmZmprnrUW3l5uboEOT09nacgewl/qz7ObrdDURSYzWaEhITAbDZDSqmuICPyVUIITJgwARaLBWfPnsWePXt0771IKZGTk4PCwkKEhoZyCbIXMVx8nOfTXnBwMIKCgtRVLb7yKZDocuLj4zFo0CA0NTVhzZo1epcDAFi7di3cbjcGDx6M3r17611Oh8Vw8WFSSlRUVAAAQkNDERAQ0GqXvt6fAol+TkBAANLT0wEA33//ve5zhXa7Hd9//z0A4Oabb4a/v7+u9XRkDBcf5+mhhIaGws/PD126dGn1dSJfJoTApEmTYDabkZeXh3379ulWi2dI7PTp0wgODsakSZM4JOZFDBcf5+m5WK1WGAwGNVw450LtRf/+/TFgwAA0NjbqumpMSomvv/4aLpcLN954I3flexnDxQdIKeFwOLB//37U1NSoja/lxL3VaoUQolW4tLxfdXU1srKyUFFRweEy8ilms1ldNbZu3TrdhsbKy8uxfv16AMBvfvMbBAQE6FJHZ8Fw8QFutxsvvPACxowZgxkzZqCxsRFA8yWNPcNf4eHhAJpDBmgOF8+RGk1NTZg5cybS09PxyCOPqJdEJvIFQgjccsstMJvNyM3N1WXVmJQSmzdvRn5+PqxWK2655RYOiXkZw6UNSSnVPy3l5eXhyy+/RENDA1avXo2jR48CaA4Nz6e8C3sutbW16iVkz5w5gzVr1qCpqQkbNmzw6Qs0Uec0cOBADBo0CI2NjVixYkWbh4vb7caXX34JRVEwatQoxMXFtenzd0YMF41dLDyA5gMoFy5ciFmzZqGoqKjVkNbOnTvVTZEOh0P9ZOdyudRjXsLCwgD8X8+lrq5OvabLgQMH1Mc3NDRgx44drWqoq6vD6tWrkZOTc9EDBC9VM5FWAgMDceeddwIA/v3vf6vXKGorp06dwubNm2EwGHD33XfDZDK16fN3RgwXDRUVFeH555/HO++8A4fDoX5dSolvvvkGf/7zn/Huu+9i9uzZrd7kd+/e3erN/eDBgwAAp9MJh8MBIYQaKl26dFGP3W9sbISUEocOHWr1/fbs2aP+v6IoeOONN3DPPffgtttuw/bt21s9l8vlwldffYV58+a1eYOnzkMIgdtuuw1hYWHIz89XrwDZFqSUWLlyJSorKxEbG4vx48dzSKwNML410tTUhJdeegmff/45jEYjIiMjMX36dAghoCgKli9frvY0vv32WxQVFSEmJgZNTU04dOgQgOaTZBVFQV5eHqSUqK+vV09E9gyHhYaGQgiBhoYGdW7lxIkTAJobsJQSx44dQ2NjIwICAlBUVIR//etfcLvdqKiowCeffIJRo0ap9/3mm2/wyCOPwOl0YufOnYiOjtbht0edQZ8+fTBu3DisWLECS5YswdSpU9tkUt3hcGD58uUAgMmTJyMqKsrrz0nsuWjGYDBg3LhxiImJgdvtxhdffKFOzFdVVbVa319aWqr+va6uDmfPngUA3HDDDQCAwsJCuFwu9aJgJpMJoaGhAJovdWw0GuF0OtHQ0AC3261e6a9fv34QQuDcuXPqMNm+fftQXFysPndWVhYqKysBNPdaFi1aBKfTibCwMNx1113cVEZeYzQa8fvf/x4mkwk7d+7E3r17vd57kVIiOzsbBw4cQGBgIO6++272WtoIw0UjBoMB06ZNw9///ncYjUYcOHAABQUFAICzZ8+itLRU3QSpKApycnIgpUR5eTkqKirg5+ennnNUVlYGh8OBmpoauFwu+Pv7qzvzLRYL/Pz84HK5UFdXB6fTiZKSEgBASkoK/P39YbfbUVZWBikldu/eDUVREB4eDn9/fxQVFSEvLw9A80KAvXv3wmg04m9/+xumTJnCQ/zIazzH8A8aNAgOhwOLFy/2ergoiqJ+0Bs6dCiGDBnCcGkjfCfRWFJSEiIjI1FVVaXOfRw/fhwNDQ3o1q0b0tLSAAA//vgjpJQoLCyEw+FAcHAwhgwZAqPRiJqaGvWP2+1GQEAAgoKCADSfMebv7w+Xy4Xq6mrU1taisrISBoMBQ4cORUhICBoaGlBQUAAppbpyLC0tDTExMXA6ndi/f7+6W7m6uhpdu3bFxIkT2ejI6ywWC/7whz9ACIFVq1bh1KlTXn2+M2fOYN26dRBCYNq0aQgMDPTq89H/YbhoLCoqCtddd53aHZdS4vDhw5BSomfPnkhOTgbQPE9SX1+PM2fOoKmpCWFhYejbty/8/f1RX18Pu92OyspK9URks9kMoHlDWmBgoHrsfnV1Nerq6mAymdC3b19ERERAURScPn0aDQ0N6nzMyJEj1R3JniE5T30DBgxAZGSkDr8t6myEELjjjjvQs2dPlJaW4osvvvBa70VKiS+//BKlpaWIjY3F5MmT+QGqDTFcNGYymZCYmAgA+OGHH9DY2IjDhw8DABISEnD99dfDaDSiuLgYFRUV6hBVdHQ0evTogaCgILhcLpSXl6OyshJSSlgsFnUuJDAwEEFBQVAUBZWVlSgrK4PT6URgYCCio6PVCfm8vDxUVlaiuLgYBoMB/fv3V+d0Dh06hPr6ejVkhg4dCqPR2Ka/J+q8unfvjnvvvRcA8MUXX7SaE9RSWVkZlixZAgCYOnUqunfv7pXnoYtjuHhBUlIShBDIzc1FUVGR2nsYMGAAevfuDbPZjJqaGhQWFuLkyZMAgF69eqFLly6wWCxwu90oKSlBeXk5gP87tBJAq/mXyspKlJSUoKmpCSEhIbBareoR4qdOnUJhYSFqampgNpvRu3dv3HDDDRBC4PTp08jPz8fJkychhMCQIUPa+DdEnZkQAr///e8RGRmJ06dPY/ny5Zr3XqSUyMjIwIkTJxAWFob7779f0+9PP4/hojEhBK6//noEBwejrKwMOTk5OHfunNp7iIyMRHh4OBobG3Hs2DGcPn0aAGCz2RAYGKjuZykpKVEPrezSpYvaszCZTOqy5PLycnVDZpcuXRAUFKTuPD579qy6JDk8PByRkZHo168f/P39UVZWhuzsbJSXlyMwMBD9+/fncAG1qfj4eNx5552QUmLRokXqa10r5eXl+PDDDyGlxB133IHrrruOr/E2xnDxgpiYGMTExMDlcmH16tWw2+0wm82Ii4uDxWJBbGwsAGDv3r0oLCwEAMTFxcHPzw8REREAmjdkejY1RkREqKu4DAaDulu/rKwM586dAwBERkbC399fDZeioiLs27cPUkrExMQgJCQEMTExCA8PR319PVauXImGhgZERUVxbwu1OSEEHn30UVitVhw/fhxfffWVZr0XKSU+//xzHD16FGFhYZgxYwZXQeqAv3EvCA4OxqBBgwAAK1euRGNjI6KiotC9e3eYTCbYbDYAwM6dO1FeXg6TyYQ+ffrAYDCga9euAJrDwTMs5gkcoLlRejaBlZaWqsudu3fvDqPRiJ49e8LPzw+VlZXYsWMHgOZekZ+fH8LCwtCnTx9IKbFu3ToAzZ8gPXtoiNqKEAIDBw7ElClToCgK3n//fZSWlmryvc+ePYsPPvgAUkrcd999GDBgAHstOmC4eIFnPT8A9WywhIQEda6kf//+AJpXbdXU1CA4OBgxMTEAgG7dugFo3kjpCZcLV3J5AqhluPTo0UP9b3BwMGpra9Wl0P369QMA+Pn5YfDgwQCaj5YBOJlP+jEYDHjyyScRERGB3NxcfPLJJxc9++5quN1uvPfeezhz5gyio6Px1FNP8fWtE4aLFwghMGbMGISEhKhfGzVqFIxGI4QQGDBgAAwGA1wuF6SUiIyMVHsnnhUthYWF6ubIrl27qp+8hBBquBQXF6vDap5wiYiIQGRkpHrwpef5PI/3HP0CNDfukSNHevvXQXRRQgj0798f999/P6SU+Mc//oEjR4784uExKSV27dqFxYsXQwiBP/7xj+jTp4+2RdMVY7h4SUJCAsaNGwegeUK+5fUjEhISEBwcrN7XZrMhODgYQgi153Lu3DmUlZW1GgbziIqKghACxcXFOHfuHIQQiI6OhhCi1aQ+AAQFBSEhIQFAc2MeNWqUGk7R0dEYNmwYhwxINwaDATNmzEDfvn1RUlKCl19+GfX19b/oe1VXV+Oll15CdXU1hgwZggcffJCvbR0xXLzEz88Pc+fOxUsvvYSPP/5YHY4Cmt/Ue/Xqpf49MTFRnXDs1q0bjEajuvveZDL9ZFgsKioKBoMBVVVVqKqqgtFoVHsuBoMBN910k3rf2NhYdcgNAHr27IlZs2Zh6NCheOmll7j2n3QXHR2Nl19+GQEBAfj3v/+NhQsXqtcqulJNTU2YO3cusrKyEBwcjDlz5qgLX0gfDBcvEUIgNjYWs2bNwq233tpqtUpwcDDS09MBNO+4b3n0SlRUlLqnBQACAgJaTeh77tPyNNnAwEC1NyKEwIQJE9Tbx40bp871AM3h8+ijj2Ljxo3qMRxEevLs2v/DH/4ARVHw1ltvYcWKFVc8/+J2u/HPf/4TCxYsgJQSjz/+OMaNG8fXts545L4XXe7F/cwzz0BRFMTFxWHEiBHq18PCwhAcHKwep2+xWNS9Lx7h4eGwWCzqNWNCQ0NbfUobOXIkZs2ahePHj+O55577SR1CCF4/nHyKn58fZs+ejePHj2PTpk149tlnYbFYcPPNN1+2HbndbixZsgQvvvgiGhoacOutt+KFF17gxcB8AP8FdOCZlH/zzTfVv3uEhobCarWqK8U8QdJSSEgIIiMjW034t1xOHBAQgOeffx5SSggh+AmOfJ4QAuHh4fjwww9x3333Yc+ePXjkkUewYMEC3H777T9Z8SWlhMPhwMKFC/HWW2/B4XAgNTUVCxcubLWQhvTDYTGdeN70L3zjN5vNrTY1xsbG/qSXERgY2GrSPj4+/ifXYRFCwGAwMFio3RBCoHfv3li8eDEGDx6MsrIyPPLII3jhhRdw6NAhdXWlw+HApk2bcPfdd2POnDmor6/HpEmT8Omnn6Jbt258zfsIhouPMZlM6unFANRlyy0ZDAYkJSWpf/ecZUbU3gkhkJCQgKVLlyItLU3tnaSlpeG2227D448/jokTJ+L222/Hhg0b4Ofnh0cffRSLFy9Gjx492A58CIfFfNCYMWPw6aefwmg0IjU19Se3CyFw1113qZduvf3229moqMMQQiAuLg7Lli3Dp59+ik8//RS5ubnYtGmTep/AwECkpqbiueeew/jx41stgiHfwHDRmBbnI916663485//jJCQEIwfP/6i3zc+Ph7/+c9/AABWq9XrV/QjulJavRZDQ0PxzDPPYPr06dixYwe+//57lJaWwmazYdKkSbjxxhvVIWO+/n0Pw0UDQgjs378fr732mqbft6amBn/72980/Z4/Z9++fewF0VXzVhtoKSQkRF3csn79eqxfv94rz8M2oA0hGfnXrLGxESdPnrzqjV++yGAwwGaz/WSBANHlsA3QhRguRESkOa4WayeklFAUhWPL1KmxHbQfDJd2Yu/evTCbzdi7d6/epRDpZu/evQgKCmI7aAcYLkREpDmGCxERaY7hQkREmmO4EBGR5hguRESkOYYLERFpjuFCRESaY7gQEZHmGC5ERKQ5hgsREWmO4UJERJpjuBARkeYYLkREpDmGCxERaY7hQkREmmO4tANSSlRWVgIAKisreaEk6pQ87aDlf8l3MVx8mN1ux/z585GQkIAJEyagsbEREyZMQEJCAubPnw+73a53iURex3bQPgnJ+PdJ69atw9SpU+FwOACg1ac0IQQAICgoCBkZGUhPT9elRiJvYztovxguPmjdunWYPHmyer3wSzEYDBBCYM2aNWxY1OGwHbRvDBcfY7fbERsbi/r6+ss2KA+DwQCz2Yz8/HxYrVbvF0jUBtgO2j/OufiYxYsXw+FwXFGDAgBFUeBwOPDZZ595uTKitsN20P6x5+JDpJRISEhAXl7eVa2EEULAZrPh+PHj6jg0UXvFdtAxMFx8SFlZGaKioq7p8RERERpWRNT22A46Bg6L+ZDa2tprenxNTY1GlRDph+2gY2C4+BCLxXJNjw8JCdGoEiL9sB10DAwXHxIREYH4+PirHi8WQiA+Ph7h4eFeqoyo7bAddAwMFx8ihMBTTz31ix779NNPcxKTOgS2g46BE/o+huv7idgOOgL2XHyM1WpFRkYGhBAwGC7/z+PZmbxixQo2KOpQ2A7aP4aLD0pPT8eaNWtgNpshhPhJN9/zNbPZjLVr12LSpEk6VUrkPWwH7RvDxUelp6cjPz8f8+bNg81ma3WbzWbDvHnzUFBQwAZFHRrbQfvFOZd2QEqJiooK1NTUICQkBOHh4Zy0pE6H7aB9YbgQEZHmOCxGRESaY7gQEZHmGC5ERKQ5hgsREWmO4UJERJpjuBARkeYYLkREpDmGCxERaY7hQkREmmO4EBGR5hguRESkOYYLERFpjuFCRESaY7gQEZHm/j8h+UhYj3axHwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 500x400 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "733a2a41",
   "metadata": {},
   "source": [
    "suggest_symbolic does not return anything that matches with it, since Bessel function isn't included in the default SYMBOLIC_LIB. We want to add Bessel to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "031db28f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  function  fitting r2   r2 loss  complexity  complexity loss  total loss\n",
      "0        0    0.000000  0.000014           0                0    0.000003\n",
      "1        x    0.001602 -0.002298           1                1    0.799540\n",
      "2      sin    0.161428 -0.253977           2                2    1.549205\n",
      "3      cos    0.161428 -0.253977           2                2    1.549205\n",
      "4    1/x^2    0.099456 -0.151116           2                2    1.569777\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "('0',\n",
       " (<function kan.utils.<lambda>(x)>,\n",
       "  <function kan.utils.<lambda>(x)>,\n",
       "  0,\n",
       "  <function kan.utils.<lambda>(x, y_th)>),\n",
       " 0.0,\n",
       " 0)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.suggest_symbolic(0,0,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4b8549a2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['x', 'x^2', 'x^3', 'x^4', 'x^5', '1/x', '1/x^2', '1/x^3', '1/x^4', '1/x^5', 'sqrt', 'x^0.5', 'x^1.5', '1/sqrt(x)', '1/x^0.5', 'exp', 'log', 'abs', 'sin', 'cos', 'tan', 'tanh', 'sgn', 'arcsin', 'arccos', 'arctan', 'arctanh', '0', 'gaussian'])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SYMBOLIC_LIB.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5db9e7cf",
   "metadata": {},
   "source": [
    "add bessel function J0 to the symbolic library. we should include a name and a pytorch implementation. c is the complexity assigned to J0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cbde1924",
   "metadata": {},
   "outputs": [],
   "source": [
    "add_symbolic('J0', torch.special.bessel_j0, c=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bda24c6d",
   "metadata": {},
   "source": [
    "After adding Bessel, we check suggest_symbolic again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "83e5cfdd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  function  fitting r2   r2 loss  complexity  complexity loss  total loss\n",
      "0        0    0.000000  0.000014           0                0    0.000003\n",
      "1       J0    0.198505 -0.319216           1                1    0.736157\n",
      "2        x    0.001602 -0.002298           1                1    0.799540\n",
      "3      sin    0.161428 -0.253977           2                2    1.549205\n",
      "4      cos    0.161428 -0.253977           2                2    1.549205\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "('0',\n",
       " (<function kan.utils.<lambda>(x)>,\n",
       "  <function kan.utils.<lambda>(x)>,\n",
       "  0,\n",
       "  <function kan.utils.<lambda>(x, y_th)>),\n",
       " 0.0,\n",
       " 0)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# J0 fitting is not very good\n",
    "model.suggest_symbolic(0,0,0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4180de14",
   "metadata": {},
   "source": [
    "The fitting r2 is still not high, this is because the ground truth is J0(20x) which involves 20 which is too large. our default search is in (-10,10). so we need to set the search range bigger in order to include 20. now J0 appears at the top of the list\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e78f4674",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  function  fitting r2   r2 loss  complexity  complexity loss  total loss\n",
      "0       J0    0.998912 -9.830484           1                1   -1.166097\n",
      "1        0    0.000000  0.000014           0                0    0.000003\n",
      "2        x    0.001602 -0.002298           1                1    0.799540\n",
      "3      cos    0.583964 -1.265186           2                2    1.346963\n",
      "4      sin    0.583964 -1.265186           2                2    1.346963\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "('J0',\n",
       " (<function torch._C._special.special_bessel_j0>,\n",
       "  J0,\n",
       "  1,\n",
       "  <function torch._C._special.special_bessel_j0>),\n",
       " 0.9989116787910461,\n",
       " 1)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.suggest_symbolic(0,0,0,a_range=(-40,40))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
